library(dplyr)
library(tidyverse)
library(brglm2)
library(brglm)
library(reshape2)
library(magrittr)
library(ggeffects)
library(geomtextpath)
library(fect)
library(readxl)
library(panelView)
library(lme4)
library(nlme)
library(modelsummary)
library(here)
library(ordinal)
library(zoo)
library(sandwich)  
library(lmtest)    
library(wesanderson)
library(patchwork)
#
# creating the autocratization onset by election ####
ert$aut_ep_onset <- ifelse(ert$year == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset[is.na(ert$aut_ep_onset)] = 0
table(ert$aut_ep_onset)
# lead 1
ert$year.lead1 <- lead(ert$year, n =1L) 
ert$aut_ep_onset1 <- ifelse(ert$year.lead1 == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset1[is.na(ert$aut_ep_onset1)] = 0
# lead 2
ert$year.lead2 <- lead(ert$year, n =2L) 
ert$aut_ep_onset2 <- ifelse(ert$year.lead2 == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset2[is.na(ert$aut_ep_onset2)] = 0
# lead 3
ert$year.lead3 <- lead(ert$year, n =3L) 
ert$aut_ep_onset3 <- ifelse(ert$year.lead3 == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset3[is.na(ert$aut_ep_onset3)] = 0
# lead 4
ert$year.lead4 <- lead(ert$year, n =4L) 
ert$aut_ep_onset4 <- ifelse(ert$year.lead4 == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset4[is.na(ert$aut_ep_onset4)] = 0
# lag 1
ert$year.lag1 <- lag(ert$year, n =1L) 
ert$aut_ep_onset.lag1 <- ifelse(ert$year.lag1 == ert$aut_ep_start_year, 1, 0)
ert$aut_ep_onset.lag1[is.na(ert$aut_ep_onset.lag1)] = 0

# TOGETHER ##
ert$aut_onset_ele <- ert$aut_ep_onset + ert$aut_ep_onset1 + 
  ert$aut_ep_onset2 + ert$aut_ep_onset3 + ert$aut_ep_onset4 + ert$aut_ep_onset.lag1
table(ert$aut_onset_ele)

ert$aut_onset_3yrs <- ert$aut_ep_onset + 
  ert$aut_ep_onset1 + ert$aut_ep_onset.lag1

ert$aut_onset_2yrs <- ert$aut_ep_onset +ert$aut_ep_onset1 

#
# creating the democratization onset by election ####
ert$dem_ep_onset <- ifelse(ert$year == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset[is.na(ert$dem_ep_onset)] = 0
# lead 1 
ert$dem_ep_onset1 <- ifelse(ert$year.lead1 == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset1[is.na(ert$dem_ep_onset1)] = 0
# lead 2 
ert$dem_ep_onset2 <- ifelse(ert$year.lead2 == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset2[is.na(ert$dem_ep_onset2)] = 0
# lead 3 
ert$dem_ep_onset3 <- ifelse(ert$year.lead3 == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset3[is.na(ert$dem_ep_onset3)] = 0
# lead 4 
ert$dem_ep_onset4 <- ifelse(ert$year.lead4 == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset4[is.na(ert$dem_ep_onset4)] = 0
# lag 1
ert$dem_ep_onset.lag1 <- ifelse(ert$year.lag1 == ert$dem_ep_start_year, 1, 0)
ert$dem_ep_onset.lag1[is.na(ert$dem_ep_onset.lag1)] = 0

# TOGETHER ##
ert$dem_onset_ele <- ert$dem_ep_onset + ert$dem_ep_onset1 + 
  ert$dem_ep_onset2 + ert$dem_ep_onset3 + ert$dem_ep_onset4+ert$dem_ep_onset.lag1
table(ert$dem_onset_ele)

ert$dem_onset_3yrs <- ert$dem_ep_onset + ert$dem_ep_onset1 + 
  ert$dem_ep_onset.lag1

ert$dem_onset_2yrs <- ert$dem_ep_onset + ert$dem_ep_onset1 

# covariates ####
vdem1 <- vdem %>% 
  group_by(e_regionpol_6C, year) %>% 
  mutate(avg.dem.reg = mean(v2x_polyarchy),
         gdp_growth = ((e_gdppc / lag(e_gdppc) - 1) * 100),
         v2xps_party_change = ((v2xps_party/lag(v2xps_party) -1)*100)
  )

vdem1$lag.v2x_polyarchy <- lag(vdem1$v2x_polyarchy, n = 1L)
vdem1$lead.v2x_polyarchy <- lead(vdem1$v2x_polyarchy, n = 1L)
vdem1$lag.avg.dem.reg <- lag(vdem1$avg.dem.reg, n = 1L)
vdem1$lag.v2xlg_legcon <- lag(vdem1$v2xlg_legcon, n = 1L)
vdem1$lag.v2elaccept <- lag(vdem1$v2elaccept, n = 1L)
vdem1$lag.v2elmulpar <- lag(vdem1$v2elmulpar, n = 1L)
vdem1$lag.e_gdppc <- lag(vdem1$e_gdppc, n = 1L)
vdem1$lag.v2pepwrsoc <- lag(vdem1$v2pepwrsoc, n = 1L)
vdem1$lag.v2xps_party <- lag(vdem1$v2xps_party, n = 1L)
vdem1$lag.v2xps_party <- lag(vdem1$v2xps_party, n = 1L)
vdem1$lag.v2x_jucon <- lag(vdem1$v2x_jucon, n = 1L)
vdem1$lag.v2elboycot <- lag(vdem1$v2elboycot, n = 1L)
vdem1$lag.v2pepwrsoc <- lag(vdem1$v2pepwrsoc, n = 1L)
vdem1$lag.v2x_diagacc <- lag(vdem1$v2x_diagacc, n = 1L)
vdem1$lag.v2x_accountability <- lag(vdem1$v2x_accountability, n = 1L)
vdem1$lag.v2smgovdom <- lag(vdem1$v2smgovdom, n = 1L)
vdem1$lag.v2x_cspart <- lag(vdem1$v2x_cspart, n = 1L)

vdem1$decade <- NA
vdem1$decade[vdem1$year >1969 & vdem1$year< 1980] <- 0
vdem1$decade[vdem1$year >1979 & vdem1$year< 1990] <- 1
vdem1$decade[vdem1$year >1989 & vdem1$year< 2000] <- 2
vdem1$decade[vdem1$year >1999 & vdem1$year< 2010] <- 3
vdem1$decade[vdem1$year >2009 & vdem1$year< 2020] <- 4
vdem1$decade <- as.factor(vdem1$decade)

## Calling the dataset created at line 168 in "psdi_computation" ##
psdi.data.clean

parsys_change <- psdi.data.clean%>%
  group_by(country_name)%>%
  mutate(psdi_change1 = ((v2ps_democracy_adj / lag(v2ps_democracy_adj) - 1) * 100),
         psdi_change2 = ((v2ps_democracy_adj / lag(v2ps_democracy_adj, n=2) - 1) * 100),
         psdi_change3 = ((v2ps_democracy_adj / lag(v2ps_democracy_adj, n=3) - 1) * 100),
         psdi_gov_change = ((v2ps_democracy_government / lag(v2ps_democracy_government) - 1) * 100),
         psdi_oppo_change = ((v2ps_democracy_opposition_adj / lag(v2ps_democracy_opposition_adj) - 1) * 100),
         ele_count = row_number()) 

parsys_change$psdi_change1[is.infinite(parsys_change$psdi_change1)] <- NA
parsys_change$psdi_change2[is.infinite(parsys_change$psdi_change2)] <- NA
parsys_change$psdi_change3[is.infinite(parsys_change$psdi_change3)] <- NA
parsys_change$psdi_gov_change[is.infinite(parsys_change$psdi_gov_change)] <- NA
parsys_change$psdi_oppo_change[is.infinite(parsys_change$psdi_oppo_change)] <- NA

## lead psdi change
parsys_change$psdi_change1_lead1 <- lead(parsys_change$psdi_change1, n =1)
parsys_change$psdi_change1_lead2 <- lead(parsys_change$psdi_change1, n =2)
parsys_change$psdi_change1_lead3 <- lead(parsys_change$psdi_change1, n =3)

## cumulative chage
parsys_change$psdi_change_cumulative2 <- parsys_change$psdi_change1 + parsys_change$psdi_change2
parsys_change$psdi_change_cumulative3 <- parsys_change$psdi_change1 + parsys_change$psdi_change2 +parsys_change$psdi_change3


parsys_change$v2ps_democracy_government <-  1 -parsys_change$v2ps_democracy_government 
parsys_change$v2ps_democracy_opposition_adj <- 1-parsys_change$v2ps_democracy_opposition_adj 

## Merges ####
# V-Dem + parsys
parsys1_vdem <- left_join(parsys_change,vdem1, by = c("country_name", "year"))
# + ERT
psdi_ert_vdem <- left_join(parsys1_vdem,ert, by = c("country_name", "year"))
# + Maddison
mpd2020 <- read_excel("mpd2020.xlsx")
mpd2020$ln_gdppc <- log(mpd2020$gdppc)
mpd2020 <- mpd2020%>%
  group_by(country_name)%>%
  mutate(gdppc_change = (gdppc / lag(gdppc) - 1) * 100)

psdi_ert_vdem_maddison <- left_join(psdi_ert_vdem,mpd2020, by = c("country_name", "year"))

# power by social group
ve_analysis <- readRDS("ve_analysis.rds")
ve_analysis$country_name <- ve_analysis$countryname
psdi_ert_vdem_maddison_marquardt <- left_join(psdi_ert_vdem_maddison,ve_analysis, by = c("country_name", "year"))

psdi_ert_vdem_maddison_marquardt$psdi_neg <- NA
psdi_ert_vdem_maddison_marquardt$psdi_neg[psdi_ert_vdem_maddison_marquardt$psdi_change1 < 0] <-  1
psdi_ert_vdem_maddison_marquardt$psdi_neg[psdi_ert_vdem_maddison_marquardt$psdi_change1 >=0] <-  0
psdi_ert_vdem_maddison_marquardt$psdi_pos <- NA
psdi_ert_vdem_maddison_marquardt$psdi_pos[psdi_ert_vdem_maddison_marquardt$psdi_change1 > 0] <-  1
psdi_ert_vdem_maddison_marquardt$psdi_pos[psdi_ert_vdem_maddison_marquardt$psdi_change1 <= 0] <-  0

## vdem_multi dataset ####
vdem_multi <- left_join(vdem1, ert, by = c("country_name", "year"))%>%
  left_join(., parsys_change, by = c("country_name", "year"))%>%
  left_join(.,mpd2020, by= c("country_name", "year"))%>%
  left_join(.,ve_analysis,by= c("country_name", "year"))

vdem_multi <- vdem_multi%>%
  filter(year >1969)

vdem_multi$trans2 <- 1 # Nothing 
vdem_multi$trans2[vdem_multi$aut_ep_onset ==1] <- 0 # Autocratization
vdem_multi$trans2[vdem_multi$dem_ep_onset ==1] <- 2 # Democratization
vdem_multi$trans2 <- as.factor(vdem_multi$trans2)
table(vdem_multi$trans2)

vdem_multi$v2ps_democracy_2 <- na.locf(vdem_multi$v2ps_democracy_adj,fromLast = TRUE,
                                       na.rm = FALSE)

# Creating variable for age of the regime
vdem_multi$reg_duration <- vdem_multi$reg_end_year-vdem_multi$reg_start_year
vdem_multi <- vdem_multi%>%
  group_by(country_name, reg_duration)%>%
  mutate(reg_count = row_number())

vdem_multi$gdppc_change[!is.finite(vdem_multi$gdppc_change)] <- NA
vdem_multi$ln_gdppc[!is.finite(vdem_multi$ln_gdppc)] <- NA

vdem_multi$ele_count_2 <- na.locf(vdem_multi$ele_count,fromLast = TRUE,
                                  na.rm = FALSE)
### ERT OUTCOME AUTOCRATIZATION ###
vdem_multi$aut_ep_outcome <- as.factor(vdem_multi$aut_ep_outcome)

vdem_multi$aut_ep_outcome_rescale <- NA
# No autocratization
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==0] <- 0
# Autocratization/Democratic Breakdown
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==1] <- 2
# Prempted Democratic Breakdown
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==2] <-3 
# Diminished Democracy
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==3] <- 4
# Averted regression
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==4] <- 5
# Regressed autocracy
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==5] <- 1
# Uncertain
vdem_multi$aut_ep_outcome_rescale[vdem_multi$aut_ep_outcome ==6] <- 6

vdem_multi$aut_ep_outcome_rescale <- as.factor(vdem_multi$aut_ep_outcome_rescale)
## Empirics ####
### REPLICATING DDLA 2013 ####
DDLA <- read_excel("DDLA.XLSX")
DDLA$country_name <- DDLA$country
DDLA$trans1 <- ifelse(DDLA$trans2 ==0, 0, 1)
parsys_ddla <-left_join(DDLA, psdi_ert_vdem_maddison_marquardt, by = c("country_name", "year"))

parsys_ddla$v2ps_democracy_adj_2 <- na.locf(parsys_ddla$v2ps_democracy_adj, fromLast = TRUE)
parsys_ddla$v2ps_democracy_opposition_adj_2 <- na.locf(parsys_ddla$v2ps_democracy_opposition_adj, fromLast = T)
parsys_ddla$v2ps_democracy_government_2 <- na.locf(parsys_ddla$v2ps_democracy_government, fromLast = T)
parsys_ddla$ele_count0 <- parsys_ddla$ele_count -1
parsys_ddla$year0 <- parsys_ddla$year - 1970

## CH 4 WITH PSDI ###
parsys_ddla$trans2 <- as.factor(parsys_ddla$trans2)
parsys_ddla$trans1 <- as.factor(parsys_ddla$trans1)
parsys_ddla$ele_count_2 <- na.locf(parsys_ddla$ele_count,fromLast = TRUE)
# Creating variable for age of the regime
parsys_ddla$reg_duration <- parsys_ddla$reg_end_year-parsys_ddla$reg_start_year
parsys_ddla <- parsys_ddla%>%
  group_by(country_name, reg_duration)%>%
  mutate(reg_count = row_number())

fmm1_psdi <- clmm(trans2 ~ gdp_1_ + ghist10 + oilmin+ indust_ +
                    ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
                    (1|year), data = parsys_ddla,
                  link = "logit",
                  na.action=na.exclude) 

fmm2_psdi <- clmm(trans2 ~ rad_all +region  +us_t  +d_nonla+
                    gdp_1_ + ghist10 + oilmin+ indust_ +
                    ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
                    (1|year), data = parsys_ddla,
                  link = "logit",
                  na.action=na.exclude) 

fmm3_psdi <- clmm(trans2 ~ rad_gov + rad_opp + region  +us_t  +d_nonla+
                    gdp_1_ + ghist10 + oilmin+ indust_ +
                    ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
                    (1|year), data = parsys_ddla,
                  link = "logit",
                  na.action=na.exclude) 

fmm4_psdi <- clmm(trans2 ~rad_all +v2ps_democracy_adj_2+
                    +region  +us_t  +d_nonla+
                    gdp_1_ + ghist10 + oilmin+ indust_ +
                    ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
                    (1|year), data = parsys_ddla,
                  link = "logit",
                  na.action=na.exclude)  

fmm5_psdi <- clmm(trans2 ~ rad_gov + rad_opp + 
                    v2ps_democracy_adj_2+
                    region  +us_t  +d_nonla+
                    gdp_1_ + ghist10 + oilmin+ indust_ +
                    ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
                    (1|year), data = parsys_ddla,
                  link = "logit",
                  na.action=na.exclude)

## Original CH 4 ###
ddla_2006 <- DDLA%>%
  filter(year <2006)
ddla_2006$trans2 <- as.factor(ddla_2006$trans2)

fmm5 <- clmm(trans2 ~ npr_all+rad_gov +rad_opp+region  +us_t  +d_nonla+
               gdp_1_ + ghist10 + oilmin+ indust_ +
               ele_count_2+ I(ele_count_2^2)+I(ele_count_2^3)+
               (1|cc_cow), data = parsys_ddla,
             link = "probit",
             na.action=na.exclude) 
summary(fmm5)

## MAIN TABLE TEXT

fmm_psdi <- list(fmm2_psdi, fmm3_psdi, fmm4_psdi, fmm5_psdi,fmm5)
msummary(fmm_psdi,estimate = "{estimate}{stars}")
#
### Figure 6. PSDI LEVELS AND REGIME CHANGE ####

psdi_aut <- glm(aut_ep_onset ~ v2ps_democracy_2+I(v2ps_democracy_2^2)+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                vdem_multi)

psdi_aut2 <- glm(aut_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 vdem_multi)

psdi_aut3 <- glm(aut_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                   lag.v2x_cspart+z+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 vdem_multi)

psdi_dem <- glm(dem_ep_onset ~ v2ps_democracy_2+I(v2ps_democracy_2^2)+
                  reg_count+ I(reg_count^2)+I(reg_count^3)+
                  (1|year), 
                family = binomial(link = "probit"), 
                vdem_multi)

psdi_dem2 <- glm(dem_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 vdem_multi)

psdi_dem3 <- glm(dem_ep_onset ~ v2ps_democracy_2 + I(v2ps_democracy_2^2)+
                   gdp_growth+gdppc_change+ln_gdppc+
                   lag.avg.dem.reg+lag.v2elaccept+lag.v2xlg_legcon+
                   lag.v2x_cspart+z+
                   reg_count+ I(reg_count^2)+I(reg_count^3)+
                   (1|year), 
                 family = binomial(link = "probit"), 
                 vdem_multi)

list_psdi_changes <- list( psdi_aut,psdi_aut2,psdi_aut3,
                           psdi_dem,psdi_dem2,psdi_dem3)
msummary(list_psdi_changes,estimate = "{estimate}{stars}", vcov = ~country_name)


## Figure 6. Predicted Probabilities ###

## CALCULATING THE VERTEX (= WHERE AUTOCRATIZATION ONSET IS MOST LIKELY) 
## ON THE RELATIONSHIP PSDI --> AUT_ONSET
aut1_psdi0 <- glm(aut_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2), 
                  family = binomial(link = "probit"), psdi_ert_vdem_maddison_marquardt)

maximum_x <- -coef(aut1_psdi0)[2] / (2 * coef(aut1_psdi0)[3])
maximum_x
maximum_y <- predict(aut1_psdi0, newdata = data.frame(v2ps_democracy_adj = maximum_x), 
                     type = "response")

aut_plot <- plot(ggeffects::ggpredict(aut1_psdi0, "v2ps_democracy_adj"), colors ="firebrick") +
  ggplot2::annotate(geom = "point", x = maximum_x, y = maximum_y, size = 3)+
  geom_textvline(label = "0.41", xintercept = 0.41, vjust = 0) +
  xlab("PSDI")+
  ylab("Pr(Autocratization Onset) %") +
  labs(title = "A")

## ON THE RELATIONSHIP PSDI --> DEM_ONSET ###
dem1_psdi0 <- glm(dem_onset_ele ~ v2ps_democracy_adj+I(v2ps_democracy_adj^2), 
                  family = binomial(link = "probit"), 
                  psdi_ert_vdem_maddison_marquardt)

maximum_dem <- -coef(dem1_psdi0)[2] / (2 * coef(dem1_psdi0)[3])
maximum_dem
maximum_demonset <- predict(dem1_psdi0, newdata = data.frame(v2ps_democracy_adj = maximum_dem), 
                            type = "response")

dem_plot <- plot(ggeffects::ggpredict(dem1_psdi0, "v2ps_democracy_adj"), colors ="deepskyblue3") +
  ggplot2::annotate(geom = "point", x = maximum_dem, y = maximum_demonset, size = 3)+
  geom_textvline(label = "0.3", xintercept = 0.3048, vjust = 0) +
  xlab("PSDI")+
  ylab("Pr(Democratization Onset) %") +
  labs(title = "B")

tiff('PSDI_aut_dem_prpro.tiff', units="in", width=10, height=5, res=600, compression = 'lzw')
aut_plot + dem_plot
dev.off()

#
### Table 2. PSDI CHANGE AUT & DEM ####

psdi_ert_vdem_maddison_marquardt$reg_duration <- 
  psdi_ert_vdem_maddison_marquardt$reg_end_year-psdi_ert_vdem_maddison_marquardt$reg_start_year

psdi_ert_vdem_maddison_marquardt_a <- psdi_ert_vdem_maddison_marquardt%>%
  group_by(country_name, reg_duration)%>%
  mutate(reg_count = row_number())

aut3_psdi_change1 <- glm(aut_onset_ele ~ psdi_change1+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)
aut3_psdi_change1_clustered <- coeftest(aut3_psdi_change1,
                                        vcov = vcovCL,
                                        type = "HC3",
                                        cluster = ~country_name)

aut3_psdi_change2 <- glm(aut_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)
aut3_psdi_change2_clustered <- coeftest(aut3_psdi_change2,
                                        vcov = vcovCL,
                                        type = "HC3",
                                        cluster = ~country_name)

aut3_psdi_change3 <- glm(aut_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)
aut3_psdi_change3_clustered <- coeftest(dem3_psdi_change3,
                                        vcov = vcovCL,
                                        type = "HC3",
                                        cluster = ~country_name)

list_psdi_changes_aut <- list(aut3_psdi_change1,aut3_psdi_change2,aut3_psdi_change3)
msummary(list_psdi_changes_aut,estimate = "{estimate}{stars}", vcov = ~country_name)
#
### PSDI CHAGE and Democratization ###
dem3_psdi_change1 <- glm(dem_onset_ele ~ psdi_change1+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)
dem3_psdi_change1_clustered <- coeftest(dem3_psdi_change1,
                                        vcov = vcovCL,
                                        type = "HC1",
                                        cluster = ~country_name)

dem3_psdi_change2 <- glm(dem_onset_ele ~ psdi_change2+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)

dem3_psdi_change2_clustered <- coeftest(dem3_psdi_change2,
                                        vcov = vcovCL,
                                        type = "HC1",
                                        cluster = ~country_name)

dem3_psdi_change3 <- glm(dem_onset_ele ~ psdi_change3+
                           reg_count+ I(reg_count^2)+I(reg_count^3)+
                           lag.avg.dem.reg+lag.v2elaccept+
                           lag.v2xlg_legcon+ gdp_growth+ln_gdppc+gdppc_change+
                           lag.v2x_cspart+z+
                           (1|year), 
                         family = binomial(link = "probit"), 
                         psdi_ert_vdem_maddison_marquardt_a)
dem3_psdi_change3_clustered <- coeftest(dem3_psdi_change3,
                                        vcov = vcovCL,
                                        type = "HC1",
                                        cluster = ~country_name)

list_psdi_changes <- list(aut3_psdi_change1,aut3_psdi_change2,aut3_psdi_change3,
                          dem3_psdi_change1,dem3_psdi_change2,dem3_psdi_change3)
msummary(list_psdi_changes,estimate = "{estimate}{stars}", vcov = ~country_name, output = "latex")


## Figure 7 ####
## AUTOCRATIZATION ###
list_clustered_aut_psdi_change <- list(aut3_psdi_change1_clustered,aut3_psdi_change2_clustered,aut3_psdi_change3_clustered)

list.psdi_changes_aut <- list("PSDI + 1" = aut3_psdi_change1_clustered,
                              "PSDI + 2" = aut3_psdi_change2_clustered,
                              "PSDI + 3" = aut3_psdi_change3_clustered)
cm_change <- c('psdi_change1' = 'PSDI Change +1',
               'psdi_change2' = 'PSDI Change +2',
               'psdi_change3' = 'PSDI Change +3') 

changes_aut <- map_dfr(c(.85, .9, .95), function(x) {
  modelplot(list.psdi_changes_aut, coef_map = cm_change, conf_level = x, draw = FALSE) %>%
    mutate(.width = x)
})

# plot it
changes_aut_plot <- ggplot(changes_aut, aes(
  y = term, x = estimate,
  xmin = conf.low, xmax = conf.high,
  color = model)) +
  ggdist::geom_pointinterval(
    position = "dodge",
    interval_size_range = c(1, 3),
    fatten_point = .1)+
  geom_vline(xintercept = 0, color = 'black')+
  theme_light() +
  xlab("Estimate") + ylab("")+
  scale_color_manual(values = wes_palette('Darjeeling1'))+
  theme(legend.position = "none")+
  labs(colour = "Model",
       title = "Autocratization")

### DEMOCRATIZATION ###
list_clustered_dem_psdi_change <- list(dem3_psdi_change1_clustered,dem3_psdi_change2_clustered,dem3_psdi_change3_clustered)

list.psdi_changes_dem <- list("PSDI + 1" = dem3_psdi_change1_clustered,
                              "PSDI + 2" = dem3_psdi_change2_clustered,
                              "PSDI + 3" = dem3_psdi_change3_clustered)
cm_change <- c('psdi_change1' = 'PSDI Change +1',
               'psdi_change2' = 'PSDI Change +2',
               'psdi_change3' = 'PSDI Change +3') 

changes_dem <- map_dfr(c(.85, .9, .95), function(x) {
  modelplot(list.psdi_changes_dem, coef_map = cm_change, conf_level = x, draw = FALSE) %>%
    mutate(.width = x)
})


# plot it
changes_dem_plot <- ggplot(changes_dem, aes(
  y = term, x = estimate,
  xmin = conf.low, xmax = conf.high,
  color = model)) +
  ggdist::geom_pointinterval(
    position = "dodge",
    interval_size_range = c(1, 3),
    fatten_point = .1)+
  geom_vline(xintercept = 0, color = 'black')+
  theme_light() +
  xlab("Estimate") + ylab("")+
  scale_color_manual(values = wes_palette('Darjeeling1'))+
  theme(legend.position = "right")+
  labs(colour = "Model",
       title = "Democratization")

tiff('PSDI_ertchange.tiff', units="in", width=10, height=6, res=600, compression = 'lzw')
changes_aut_plot + changes_dem_plot
dev.off()


